Import des données : matrice d’expression

On a, pour chaque gène et chaque condition, son niveau d’expression en sortie de quantification. On labelle les conditions avec le code suivant : lettre majuscule pour le niveau fort, minuscule pour le niveau faible. Le réplicat est donné après l’underscore.

# quantification file
data <- read.csv("quantifFiles/quantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, grepl("R", colnames(data))]

getLabel("R6")

[1] "cnF_3"

getExactCondition("R4")

[1] At_AmbientCO2_LowNitrate_Fe1
48 Levels: At_AmbientCO2_HighNitrate_Fe1 ... Sl_ElevatedCO2_LowNitrate_FeStarvation3

getCondition("R4")

[1] "At_AmbientCO2_LowNitrate_Fe"

keep <- rowSums(data) >= 10
data <- data[keep, ]
group <- sapply(colnames(data), getLabel, with.rep = F)
colnames(data) <- sapply(colnames(data), getLabel)
head(data)

          cNF_3 cNF_2 cNF_1 cnF_2 cnF_1 cnF_3 CNF_1 CnF_2 CnF_1 CnF_3 cNf_1
AT1G01010  1526  1006  1116  1275   967  1018   854  1132  1294  1364  2325
AT1G01020   416   285   289   349   364   370   226   513   502   561   461
AT1G01030    31    15    19    29    36    28    12    47    34    47    18
          cnf_2 cnf_1 cNf_2 cNf_3 cnf_3 Cnf_3 CNf_1 Cnf_1 Cnf_2 CNf_3 CNf_2
AT1G01010  2113  2193  2564  2364  2074  1987  2027  1754  1697  1537  1898
AT1G01020   407   432   614   380   502   484   467   426   415   413   462
AT1G01030    40    32    44    37    27    42    39    36    40    37    37
          CNF_3 CNF_2
AT1G01010   816   912
AT1G01020   223   312
AT1G01030    15    19
 [ reached 'max' / getOption("max.print") -- omitted 3 rows ]

dim(data)

[1] 23342    24

getExpression("AT1G01010")

 [1] "cNF_3" "cNF_2" "cNF_1" "cnF_2" "cnF_1" "cnF_3" "CNF_1" "CnF_2" "CnF_1"
[10] "CnF_3" "cNf_1" "cnf_2" "cnf_1" "cNf_2" "cNf_3" "cnf_3" "Cnf_3" "CNf_1"
[19] "Cnf_1" "Cnf_2" "CNf_3" "CNf_2" "CNF_3" "CNF_2"
    cNF_3     cNF_2     cNF_1     cnF_2     cnF_1     cnF_3     CNF_1     CnF_2 
1404.4841 1137.1461 1370.4972 1234.5941  975.1056 1079.3725 1464.9136 1068.1962 
    CnF_1     CnF_3     cNf_1     cnf_2     cnf_1     cNf_2     cNf_3     cnf_3 
1184.6501 1186.8721 2176.3477 2035.5245 1879.4058 2023.3233 2118.1681 1941.6821 
    Cnf_3     CNf_1     Cnf_1     Cnf_2     CNf_3     CNf_2     CNF_3     CNF_2 
1766.0503 1878.6843 1666.8990 1566.9834 1652.5818 1823.3202 1342.1000 1267.1018 

Effet nitrate en conditions contrôle

On définit les conditions contrôle comme suit : CO2 ambiant et fort fer.

g = list()
method = "edger"
labels <- c("cNF", "cnF")
genes1 <- dualDE(data, labels, pval = 0.01, method = method)

  (Intercept) groupcnF
1           1        0
2           1        0
3           1        0
4           1        1
5           1        1
6           1        1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

[1] "cNF_3" "cNF_2" "cNF_1" "cnF_2" "cnF_1" "cnF_3"
    cNF_3     cNF_2     cNF_1     cnF_2     cnF_1     cnF_3 
0.9859482 0.9865531 0.9545113 1.0150637 1.0299591 1.0279646 

[1] "2803  genes DE"

g[[paste(labels[1], labels[2])]] = as.vector(genes1)

Fort CO2

labels <- c("CNF", "CnF")
genes2 <- dualDE(data, labels, pval = 0.01)

  (Intercept) groupCnF
1           1        0
2           1        0
3           1        0
4           1        1
5           1        1
6           1        1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

[1] "CNF_1" "CNF_3" "CNF_2" "CnF_2" "CnF_1" "CnF_3"
    CNF_1     CNF_3     CNF_2     CnF_2     CnF_1     CnF_3 
0.9793780 0.9377198 0.9911882 1.0424709 1.0209230 1.0283201 

[1] "1811  genes DE"

g[[paste(labels[1], labels[2])]] = genes2

Low iron

labels <- c("cNf", "cnf")
genes3 <- dualDE(data, labels, pval = 0.01)

  (Intercept) groupcnf
1           1        0
2           1        0
3           1        0
4           1        1
5           1        1
6           1        1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

[1] "cNf_1" "cNf_2" "cNf_3" "cnf_2" "cnf_1" "cnf_3"
    cNf_1     cNf_2     cNf_3     cnf_2     cnf_1     cnf_3 
0.9927341 1.0205843 0.9891514 0.9943126 1.0031887 1.0000289 

[1] "5356  genes DE"

g[[paste(labels[1], labels[2])]] = genes3

Low iron and high CO2

labels <- c("CNf", "Cnf")
genes4 <- dualDE(data, labels, pval = 0.01)

  (Intercept) groupCnf
1           1        0
2           1        0
3           1        0
4           1        1
5           1        1
6           1        1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

[1] "CNf_1" "CNf_3" "CNf_2" "Cnf_3" "Cnf_1" "Cnf_2"
    CNf_1     CNf_3     CNf_2     Cnf_3     Cnf_1     Cnf_2 
1.0043575 1.0334462 1.0275269 0.9519021 0.9851404 0.9976269 

[1] "2673  genes DE"

g[[paste(labels[1], labels[2])]] = genes4

Venn diagram

On visualise les gènes différentiellement exprimés en commun entre les différents niveaux des autres facteurs.

library(ggVennDiagram)
library(VennDiagram)

gene_list <- list()
for (comp in names(g)) {
    print(g[[comp]])
    gene_list[[comp]] <- g[[comp]]$gene_id
}

    gene_id   a.value   m.value       p.value       q.value rank estimatedDEG
1 AT1G29100  7.527717 -3.844547 3.578192e-138 8.352216e-134    1            1
2 AT1G19250  8.839658 -3.648727 1.878457e-127 2.192347e-123    2            1
3 AT2G13810  8.469376 -4.131913 1.462868e-124 1.138209e-120    3            1
4 AT5G33355  6.962703 -5.152094 8.806368e-116 5.138956e-112    4            1
5 AT4G21680  7.657504 -3.886983 1.141541e-114 5.329171e-111    5            1
6 AT4G26050 10.258910 -2.260091  1.152467e-95  4.483480e-92    6            1
7 AT2G41810  7.414076 -4.047054  8.398511e-92  2.800544e-88    7            1
8 AT4G22640  9.630374 -1.975045  1.839254e-80  5.366485e-77    8            1
9 AT4G25100  9.022566 -2.067878  7.469292e-71  1.937202e-67    9            1
  upreg
1     0
2     0
3     0
4     0
5     0
6     0
7     0
8     0
9     0
 [ reached 'max' / getOption("max.print") -- omitted 2794 rows ]
    gene_id  a.value   m.value       p.value       q.value rank estimatedDEG
1 AT2G41810 6.909960 -4.891457 5.753358e-166 1.342949e-161    1            1
2 AT4G21680 7.158203 -3.903113 1.079396e-127 1.259763e-123    2            1
3 AT5G33355 7.532822 -3.942011  1.137090e-79  8.847318e-76    3            1
4 AT4G26050 9.901752 -2.410582  3.640522e-76  2.124427e-72    4            1
5 AT2G43890 4.534486 -4.585951  5.500546e-52  2.567875e-48    5            1
6 AT4G13420 7.045056  3.418569  6.910049e-50  2.688239e-46    6            1
7 AT4G11393 2.688467 -6.030749  5.311806e-40  1.771260e-36    7            1
8 AT3G45060 3.835388 -4.139280  1.476109e-35  4.306917e-32    8            1
9 AT1G51840 8.737173  1.711245  6.237559e-35  1.617746e-31    9            1
  upreg
1     0
2     0
3     0
4     0
5     0
6     1
7     0
8     0
9     1
 [ reached 'max' / getOption("max.print") -- omitted 1802 rows ]
    gene_id   a.value   m.value       p.value       q.value rank estimatedDEG
1 AT1G08100  8.250556  7.543363 1.061212e-240 2.477082e-236    1            1
2 AT4G19680  7.581300 -5.390613 1.732236e-215 2.021692e-211    2            1
3 AT4G31970 10.071433  4.120985 6.675358e-196 5.193874e-192    3            1
4 AT4G13250 10.098952  3.612679 1.457726e-189 8.506561e-186    4            1
5 AT1G60050  9.440827  5.561353 2.298538e-189 1.073050e-185    5            1
6 AT4G01390  8.024264  5.780911 9.500326e-178 3.695943e-174    6            1
7 AT1G26390  9.011257  3.803360 1.268011e-153 4.228272e-150    7            1
8 AT1G60470  9.278972  3.398056 1.922528e-150 5.609457e-147    8            1
9 AT1G52070  6.816105 -6.878196 7.952279e-147 2.062468e-143    9            1
  upreg
1     1
2     0
3     1
4     1
5     1
6     1
7     1
8     1
9     0
 [ reached 'max' / getOption("max.print") -- omitted 5347 rows ]
    gene_id   a.value   m.value       p.value       q.value rank estimatedDEG
1 AT5G37990  8.804797  4.793179 1.242358e-182 2.899912e-178    1            1
2 AT3G55970  6.219576  5.310631  9.095523e-97  1.061538e-92    2            1
3 AT2G43510  9.062836  3.025254  2.217784e-86  1.725584e-82    3            1
4 AT4G21680 11.255880 -2.409529  8.268693e-85  4.825196e-81    4            1
5 AT1G32350  7.026600  3.559752  4.914056e-73  2.294078e-69    5            1
6 AT3G45060  9.354052 -2.056677  4.116341e-69  1.601394e-65    6            1
7 AT1G26390  9.590338  2.992967  4.952257e-68  1.651365e-64    7            1
8 AT2G29460  9.917335  1.860758  3.822791e-66  1.115395e-62    8            1
9 AT1G52890  8.325089  3.151416  1.724274e-65  4.472000e-62    9            1
  upreg
1     1
2     1
3     1
4     0
5     1
6     0
7     1
8     1
9     1
 [ reached 'max' / getOption("max.print") -- omitted 2664 rows ]

ggVennDiagram(gene_list)

partitions <- get.venn.partitions(gene_list, keep.elements = T)
partitions$shared <- rowSums(partitions[1:4])
partitions <- partitions[order(-partitions$shared), ]

common_genes <- unlist(partitions[1, "..values.."])
results <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "description", 
    "external_gene_name", "entrezgene_id"), values = common_genes, mart = mart)
results <- results[!rownames(results) %in% which(duplicated(results$ensembl_gene_id)), 
    ]
kable(results)

sharedBy4 <- unique(unlist(subset(partitions, partitions$shared == 3)$..values..))
a <- OntologyProfileAt(sharedBy4)

Summary graph

d <- data.frame(matrix(ncol = 4, nrow = 2))
colnames(d) <- names(g)
row.names(d) <- c("up", "down")
for (comp in names(g)) {
    d["up", comp] <- sum(g[[comp]]$upreg == 1)
    d["down", comp] <- sum(g[[comp]]$upreg == 0)
}
res <- melt(d)
res$reg = rep(c("up", "down"), 4)


ggplot(res, aes(fill = reg, y = value, x = variable)) + geom_bar(position = "stack", 
    stat = "identity") + ggtitle("Nitrate effet on gene regulation") + xlab("") + 
    ylab("Number of differentially expressed genes") + coord_flip()